!-------------------------------------------------------------------------------
! FileName:     gol_mpi_virtual.F90
! 
! Compile :     $ mpif90 -g -o gol_mpi_virtual.x -Wall gol_mpi_virtual.F90
!
! Run     :     $ mpirun -np 4 -hostfile hosts gol_mpi_virtual.x
!
! Be careful:   1. GridX, GridY must be same size. (ex) 10x10, 20x20, 100x100
!               2. To do preprocessor, .F90 is required.
!               3. # of processor is 4
!
! Tab size:     8
! Version :     1.09
!
! Updated by    2011. 11. 01
!-------------------------------------------------------------------------------

PROGRAM gol_mpi_virtual
IMPLICIT NONE
INCLUDE 'mpif.h'
        INTEGER :: ROOT = 0
        INTEGER nRank, nProcs, iErr, status(MPI_STATUS_SIZE), row_type, comm_cart
        INTEGER i, j, k, l, step
        INTEGER im, ip, jm, jp, nsum
        INTEGER i1, i2, i2m, j1, j2, j2m
        INTEGER ninom, njnom                    ! ig, jg
        INTEGER, DIMENSION(16) :: request
        REAL    arand
        INTEGER sqrt_nprocs, nLiveCell

        INTEGER, PARAMETER ::   nGridXSize=20, nGridYSize=20, nSteps=100, nSleep=1
        INTEGER, DIMENSION(:, :), ALLOCATABLE   :: global_buf
        INTEGER, DIMENSION(:, :), ALLOCATABLE   :: old, new

        INTEGER :: ndims = 2, reorder, mycart_rank, ew_dir, ns_dir
        INTEGER :: proc_east, proc_west, proc_north, proc_south
        INTEGER :: proc_ne, proc_nw, proc_se, proc_sw
        INTEGER dim_size(0:1), periods(0:1)
        
        INTEGER i1n, i2n, j1n, j2n
        INTEGER i1p, j1p, nitot
        INTEGER i1p_local, i2m_local, j1p_local, j2m_local

        ! MPI Start     
        CALL MPI_INIT(iErr)
        CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nProcs, iErr)
        CALL MPI_COMM_RANK(MPI_COMM_WORLD, nRank, iErr)

        dim_size(0) = nProcs / ndims    ! dim_size(0) = 4 / 2 = 2
        dim_size(1) = nProcs / ndims    ! dim_size(10 = 4 / 2 = 2

        ninom = nGridXSize / dim_size(0)! ninom = 20 / 2 = 10
        njnom = nGridYSize / dim_size(1)! njnom = 20 / 2 = 10

        ! Create virtual topology
        ns_dir = 0
        ew_dir = 1

        periods(0) = 1
        periods(1) = 1
        reorder = 0

        sqrt_nprocs = SQRT(REAL(nProcs))

        CALL MPI_CART_CREATE(MPI_COMM_WORLD, ndims, dim_size, periods, reorder, comm_cart, iErr)
        CALL MPI_CART_SHIFT(comm_cart, ew_dir, 1, mycart_rank, proc_east, iErr)
        CALL MPI_CART_SHIFT(comm_cart, ew_dir, -1, mycart_rank, proc_west, iErr)
        CALL MPI_CART_SHIFT(comm_cart, ns_dir, -1, mycart_rank, proc_north, iErr)
        CALL MPI_CART_SHIFT(comm_cart, ns_dir, 1, mycart_rank, proc_south, iErr)

        ! Virtual topology has been used to define adjacent processors
        ! to the N, S, E, W: still have to define processors on corners
        SELECT CASE (nRank)
                CASE (0)
                        proc_nw = 3
                        proc_ne = 3
                        proc_se = 3
                        proc_sw = 3
                CASE (1)
                        proc_nw = 2
                        proc_ne = 2
                        proc_se = 2
                        proc_sw = 2
                CASE (2)
                        proc_nw = 1
                        proc_ne = 1
                        proc_se = 1
                        proc_sw = 1
                CASE (3)
                        proc_nw = 0
                        proc_ne = 0
                        proc_se = 0
                        proc_sw = 0
        END SELECT

        ! Add in Fortran
        ! Norminal starting and ending indices (without ghost cells)
        i1n = (nRank / 2) * ninom + 1   ! (nRank 0/2) * 10 + 1 = 1
                                        ! (nRank 1/2) * 10 + 1 = 1
                                        ! (nRank 2/2) * 10 + 1 = 11
                                        ! (nRank 3/2) * 10 + 1 = 11
                                                
        i2n = i1n + ninom - 1           ! nRank(0) 1 + 10 - 1 = 10
                                        ! nRank(1) 1 + 10 - 1 = 10
                                        ! nRank(2) 11 + 10 - 1 = 20
                                        ! nRank(3) 11 + 10 - 1 = 20
        
        j1n = MOD(nRank, 2) * njnom + 1 ! nRank(0) 0 * 10 + 1 = 1
                                        ! nRank(1) 1 * 10 + 1 = 11
                                        ! nRank(2) 0 * 10 + 1 = 1
                                        ! nRank(3) 1 * 10 + 1 = 11
                                        
        j2n = j1n + njnom - 1           ! nRank(0) 1 + 10 - 1 = 10
                                        ! nRank(1) 11 + 10 - 1 = 20
                                        ! nRank(2) 1 + 10 - 1 = 10
                                        ! nRank(3) 11 + 10 - 1 = 20
        
        ! Local starting and ending indices, including 2 ghost cells
        i1 = i1n - 1                   ! nRank(0) 1 - 1 = 0
                                       ! nRank(1) 1 - 1 = 0
                                       ! nRank(2) 11 - 1 = 10
                                       ! nRank(3) 11 - 1 = 10
                                       
        i2 = i2n + 1                   ! nRank(0) 10 + 1 = 11
                                       ! nRank(1) 10 + 1 = 11
                                       ! nRank(2) 20 + 1 = 21
                                       ! nRank(3) 20 + 1 = 21
                                               
        i1p = i1 + 1                   ! nRank(0) 0 + 1 = 1
                                       ! nRank(1) 0 + 1 = 1
                                       ! nRank(2) 10 + 1 = 11
                                       ! nRank(3) 10 + 1 = 11
                                              
        i2m = i2 - 1                   ! nRank(0) 11 - 1 = 10
                                       ! nRank(1) 11 - 1 = 10
                                       ! nRank(2) 21 - 1 = 20
                                       ! nRank(3) 21 - 1 = 20
        
        j1 = j1n - 1                   ! nRank(0) 1 - 1 = 0
                                       ! nRank(1) 11 - 1 = 10
                                       ! nRank(2) 1 - 1 = 0
                                       ! nRank(3) 11 - 1 = 10
       
        j2 = j2n +1                    ! nRank(0) 10 + 1 = 11
                                       ! nRank(1) 20 + 1 = 21
                                       ! nRank(2) 10 + 1 = 11
                                       ! nRank(3) 20 + 1 = 21
           
        j1p = j1 + 1                   ! nRank(0) 0 + 1 = 1
                                       ! nRank(1) 10 + 1 = 11
                                       ! nRank(2) 0 + 1 = 1
                                       ! nRank(3) 10 + 1 = 11
       
        j2m = j2 - 1                   ! nRank(0) 11 - 1 = 10
                                       ! nRank(1) 21 - 1 = 20
                                       ! nRank(2) 11 - 1 = 10
                                       ! nRank(3) 21 - 1 = 20
       
        nitot = i2 - i1 + 1            ! nRank(0) 11 - 0 + 1 = 12
                                       ! nRank(1) 11 - 0 + 1 = 12
                                       ! nRank(2) 21 - 10 + 1 = 12
                                       ! nRank(3) 21 - 10 + 1 = 12

        ! Allocate Memory
        ALLOCATE(old(i1:i2, j1:j2))  ! nRank(0)      0:11, 0:11
                                     ! nRank(1)      0:11, 10:21
                                     ! nRank(2)      10:21, 0:11
                                     ! nRank(3)      10:21, 10:21
                                             
        ALLOCATE(new(i1:i2, j1:j2))
        
        IF (nRank == ROOT) THEN
                ALLOCATE(global_buf(1:nGridXSize, 1:nGridYSize))  ! nGridXSize=nGridYSize=20
        END IF

        ! Initialize elements of old to 0 or 1.  We're doing some
        ! sleight of hand here to make sure we initialize to the
        ! same values as in the serial case. The random_number
        ! function is called for every i and j, even if they are
        ! not on the current processor, to get the same random
        ! distribution as the serial case, but they are only used
        ! if this i and j reside on the current procesor.
        
        ! Get Random #               
        CALL INIT_RANDOM_SEED()
     
        DO i=1, nGridXSize
                DO j=1, nGridYSize
                        CALL RANDOM_NUMBER(arand)
                        IF (i > i1 .AND. i<i2 .AND. j>j1 .AND. j<j2) THEN
                                old(i, j) = NINT(arand)
                        END IF
                END DO
        END DO
        
        ! Create derived type for single row of array.
        ! There are nj "blocks," each containing 1 element,
        ! with a stride of nitot between the blocks
        CALL MPI_TYPE_CONTIGUOUS(njnom, MPI_INTEGER, row_type, iErr)
        CALL MPI_TYPE_COMMIT(row_type, iErr)

        ! ----------------------------------------------------------------------
        ! Main Loop : time steps
        ! ----------------------------------------------------------------------
        DO step=1, nSteps
                ! Use derived type "row_type" to transfer row
                ! send and receive
                
                ! top, bottom
                CALL MPI_ISEND(old(i1p, j1p), 1, row_type, proc_north, 0, comm_cart, request(1), iErr)
                CALL MPI_IRECV(old(i2, j1p), 1, row_type, proc_south, 0, comm_cart, request(2), iErr)
                CALL MPI_ISEND(old(i2m, j1p), 1, row_type, proc_south, 1, comm_cart, request(3), iErr)
                CALL MPI_IRECV(old(i1, j1p), 1, row_type, proc_north, 1, comm_cart, request(4), iErr)
        
                ! left, right
                CALL MPI_ISEND(old(i1p, j1p), ninom, MPI_INTEGER, proc_west, 2, comm_cart, request(5), iErr)
                CALL MPI_IRECV(old(i1p, j2), njnom, MPI_INTEGER, proc_east, 2, comm_cart, request(6), iErr)
                CALL MPI_ISEND(old(i1p, j2m), njnom, MPI_INTEGER, proc_east, 3, comm_cart, request(7), iErr)
                CALL MPI_IRECV(old(i1p, j1), njnom, MPI_INTEGER, proc_west, 3, comm_cart, request(8), iErr)

                ! 4 corner
                CALL MPI_ISEND(old(i1p, j1p), 1, MPI_INTEGER, proc_nw, 4, comm_cart, request(9), iErr)
                CALL MPI_IRECV(old(i2, i2), 1, MPI_INTEGER, proc_se, 4, comm_cart, request(10), iErr)
                CALL MPI_ISEND(old(i1p, j2m), 1, MPI_INTEGER, proc_ne, 5, comm_cart, request(11), iErr)
                CALL MPI_IRECV(old(i2, j1), 1, MPI_INTEGER, proc_sw, 5, comm_cart, request(12), iErr)

                CALL MPI_ISEND(old(i2m, j1p), 1, MPI_INTEGER, proc_sw, 6, comm_cart, request(13), iErr)
                CALL MPI_IRECV(old(i1, j2), 1, MPI_INTEGER, proc_ne, 6, comm_cart, request(14), iErr)
                CALL MPI_ISEND(old(i2m, j2m), 1, MPI_INTEGER, proc_se, 7, comm_cart, request(15), iErr)
                CALL MPI_IRECV(old(i1, j1), 1, MPI_INTEGER, proc_nw, 7, comm_cart, request(16), iErr)

                DO i=1, 16
                        CALL MPI_WAIT(request(i), status, iErr)
                END DO
                
                ! Update states of cells
                DO j=j1p, j2m
                        DO i=i1p, i2m
                                ip = i + 1
                                im = i - 1
                                jp = j + 1
                                jm = j - 1

                                nsum = old(im, jp) + old(i, jp) + old(ip, jp) + old(im, j) + old(ip, j) + old(im, jm) + old(i, jm) + old(ip, jm)
                                
                                SELECT CASE (nsum)
                                        CASE (3)
                                                new(i, j) = 1
                                        CASE (2)
                                                new(i, j) = old(i, j)
                                        CASE DEFAULT
                                                new(i, j) = 0
                                END SELECT
                        END DO
                END DO

                ! copy new state -> old state
                old(i1p:i2m, j1p:j2m) = new (i1p:i2m, j1p:j2m)
        
                !Send my own do-work data to ROOT(0)
                CALL MPI_SEND(i1p, 1, MPI_INTEGER, ROOT, 51, MPI_COMM_WORLD, iErr)
                CALL MPI_SEND(i2m, 1, MPI_INTEGER, ROOT, 52, MPI_COMM_WORLD, iErr)
                CALL MPI_SEND(j1p, 1, MPI_INTEGER, ROOT, 53, MPI_COMM_WORLD, iErr)
                CALL MPI_SEND(j2m, 1, MPI_INTEGER, ROOT, 54, MPI_COMM_WORLD, iErr)
                        
                DO j=j1p, j2m
                        CALL MPI_SEND(old(i1p, j), 1, row_type, ROOT, 55, MPI_COMM_WORLD, iErr)
                END DO
                
                IF (nRank == ROOT) THEN
                        DO k=0, nProcs-1
                                CALL MPI_RECV(i1p_local, 1, MPI_INTEGER, k, 51, MPI_COMM_WORLD, status, iErr)
                                CALL MPI_RECV(i2m_local, 1, MPI_INTEGER, k, 52, MPI_COMM_WORLD, status, iErr)
                                CALL MPI_RECV(j1p_local, 1, MPI_INTEGER, k, 53, MPI_COMM_WORLD, status, iErr)
                                CALL MPI_RECV(j2m_local, 1, MPI_INTEGER, k, 54, MPI_COMM_WORLD, status, iErr)

                                DO l=j1p_local, j2m_local
                                        CALL MPI_RECV(global_buf(i1p_local, l), 1, row_type, k, 55, MPI_COMM_WORLD, status, iErr)
                                END DO
                        END DO
                END IF
                
                CALL MPI_BARRIER(MPI_COMM_WORLD, iErr)
                
                IF (nRank == ROOT) THEN
                        CALL showGridState(global_buf, 1, 1, nGridXSize, nGridYSize, 20)
                        
                        WRITE(*, *)
                        WRITE(*, '(A, I4, A, I6)') 'msg> nStep = ', step, ' /', nSteps
                        nLiveCell = 0
                        
                        DO j=1, nGridXSize
                                DO k=1, nGridYSize
                                        nLiveCell = nLiveCell + global_buf(j, k)
                                END DO
                        END DO
                        
                        WRITE(*, '(A, I4, A, I4)') 'msg> nLiveCell = ', nLiveCell, ' /', nGridXSize*nGridYSize
                        WRITE(*, '(A, I2, A)') 'msg> Wait a ', nSleep, ' sec'
                        WRITE(*, *)
                        
                        CALL Sleep(nSleep)
                END IF
        END DO
        ! ----------------------------------------------------------------------
        ! End of step
        ! ----------------------------------------------------------------------

        ! Free memory
        DEALLOCATE(old)
        DEALLOCATE(new)
                
        IF (nRank == ROOT) THEN
                DEALLOCATE(global_buf)
        END IF

        ! MPI End
        CALL MPI_FINALIZE(iErr)

!------------------------------------------------------------
CONTAINS
SUBROUTINE showGridState(gol_buf, i1, j1, i2, j2, nCol)
INTEGER i1, i2, j1, j2, nCol
INTEGER, DIMENSION(i1:i2, j1:j2)                :: gol_buf
        INTEGER ii, jj
        
        WRITE(*, 400, advance='no')     '   '
        
        DO ii=1, nCol
                WRITE(*, 100, advance='no')     ii
        END DO
        
        WRITE(*, *)                             ! new line
        WRITE(*, 400, advance='no')     '   '
        
        Do ii=1, nCol
                WRITE(*, 200, advance='no')     '--'
        END DO
        
        WRITE(*, *)                             ! new line
                
100     FORMAT(I2, X)
200     FORMAT(A2, X)

        DO ii=i1, i2
                WRITE(*, 300, advance='no')     ii, '|'
                
                DO jj=j1, j2
                        IF (gol_buf(ii, jj) == 0) THEN
                                WRITE(*, 400, advance='no')     ' - '
                        ELSE IF (gol_buf(ii, jj) == 1) THEN
                                WRITE(*, 400, advance='no')     ' # '
                        ELSE
                                WRITE(*, 400, advance='no') ' $ '
                        END IF
                END DO
                WRITE(*, *)                             ! new line
        END DO 

300     FORMAT(I2, A)
400     FORMAT(A)
END SUBROUTINE

SUBROUTINE INIT_RANDOM_SEED()
            INTEGER :: i, n, clock
            INTEGER, DIMENSION(:), ALLOCATABLE :: seed
          
            CALL RANDOM_SEED(size = n)
            ALLOCATE(seed(n))
          
            CALL SYSTEM_CLOCK(COUNT=clock)
          
            seed = clock + 37 * (/ (i - 1, i = 1, n) /)
            CALL RANDOM_SEED(PUT = seed)
          
            DEALLOCATE(seed)
END SUBROUTINE

!-------------------------------------------------------------------------------
END
